#include <stdio.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <limits.h>
#include <values.h>
/* #include "/home/njozefow/working/solver/concorde/concorde.h" */ /* Mettre le chemin de l INSA */
#include "/gei/concorde/src/concorde.h"
#include "/usr/include/lpsolve/lp_lib.h" /* Mettre le chemin de l INSA */
#include "tsp_data.h"
#include "bacut.h"

#define INTEGER_AND_FEASIBLE -100
#define DOMINATED -200
#define EPSILON 0.000000000001

void bacut_init(TSP_DATA *d, BACUT *bc){
  bc->vect  = (double *)malloc((1+(d->nbn*(d->nbn-1)/2))*sizeof(double));
  bc->colno = (int *)malloc((1+(d->nbn*(d->nbn-1)/2))*sizeof(int));

  bc->x     = (double *)malloc((d->nbn*(d->nbn-1)/2)*sizeof(double));
  
  bc->bd    = (int *)malloc((d->nbn*(d->nbn-1)/2)*sizeof(int));
 
  bc->elist = (int *)malloc(d->nbn*(d->nbn-1)*sizeof(int));

  bc->compscount = NULL;
  bc->comps      = NULL;
  bc->ncomp      = 0;

  bc->s = (int *)malloc(d->nbn*sizeof(int));

  bc->dlen = (double *)malloc((d->nbn*(d->nbn-1)/2)*sizeof(double));
}

void bacut_free(BACUT *bc){
  free(bc->vect);
  free(bc->colno);

  free(bc->x);

  free(bc->bd);

  free(bc->elist);

  free(bc->s);

  free(bc->dlen);
}
	    
void bacut_build_lp(TSP_DATA *d, BACUT *bc){
  int i, j, ptr;

  /* Nombre de contraintes et nombre de variables */
  bc->lp = make_lp(0, d->nbn*(d->nbn-1)/2);

  /* Construction de la fonction objectif */
  ptr = 1;
  for (i = 0; i < d->nbn-1; i++)
    for (j = i+1; j < d->nbn; j++)
      bc->vect[ptr++] = d->d[i][j];
  set_obj_fn(bc->lp, bc->vect); 
  set_minim(bc->lp);

  /* Ajout des contraintes de degre */
  for (i = 0; i < d->nbn; i++){
    ptr = 0;
    
    for (j = 0; j < i; j++){
      bc->vect[ptr]    = 1.0;
      bc->colno[ptr++] = j*(2*d->nbn - 1 - j)/2 + (i - j - 1) + 1;
    }
  
    for (j = i+1; j < d->nbn; j++){
      bc->vect[ptr]    = 1.0;
      bc->colno[ptr++] = i*(2*d->nbn - 1 - i)/2 + (j - i - 1) + 1;
    }

    add_constraintex(bc->lp, ptr, bc->vect, bc->colno, EQ, 2.0);
  }

  /* Borne sur les variables */
  for (i = 0; i < d->nbn*(d->nbn-1)/2; i++)
    set_bounds(bc->lp, i+1, 0.0, 1.0);

  /* print_lp(bc->lp); */
}

void bacut_set_bounds(TSP_DATA *d, BACUT *bc){
  int i;

  for (i = 0; i < d->nbn*(d->nbn-1)/2; i++)
    if (bc->bd[i] == -1)
      set_bounds(bc->lp, i+1, 0.0, 1.0);
    else if (bc->bd[i] == 0)
      set_bounds(bc->lp, i+1, 0.0, 0.0);
    else if (bc->bd[i] == 1)
      set_bounds(bc->lp, i+1, 1.0, 1.0);
}

int bacut_is_integer_and_feasible(TSP_DATA *d, BACUT *bc){
  int i, j, ptr, ncomp;

  /* Tester si le graphe est connexe */
  ptr = 0;
  bc->ecount = 0;
  for (i = 0; i < d->nbn-1; i++)
    for (j = i+1; j < d->nbn; j++)
      if ((float)bc->x[ptr++] > EPSILON){
	bc->elist[2*(bc->ecount)]   = i;
	bc->elist[2*(bc->ecount)+1] = j;
	bc->ecount++;
      }

  CCcut_connect_components(d->nbn, bc->ecount, bc->elist, NULL, &(bc->ncomp), &(bc->compscount), &(bc->comps));

  /* Tester si la solution est entiere */
  for (i = 0; i < d->nbn*(d->nbn-1)/2; i++)
    if (((float)bc->x[i] > EPSILON) && ((float)bc->x[i] < 1.0))
      return 0;
  
  if (bc->ncomp != 1) return 0;

  return 1;
}

int bacut_solve_lp(TSP_DATA *d, BACUT *bc){
  int ret;
  
  ret = solve(bc->lp);

  if (ret == OPTIMAL){
    get_variables(bc->lp, bc->x);

    bc->lb = get_objective(bc->lp);
  
    if (bc->lb >= bc->best_length) ret = DOMINATED;

    if (bacut_is_integer_and_feasible(d, bc)){
      if ((float)bc->lb < (float)bc->best_length)
	bc->best_length = bc->lb;

      ret = INTEGER_AND_FEASIBLE;
    }
  }
  
  return ret;
}

/* La fonction a modifier */
int bacut_constraint_generation(TSP_DATA *d, BACUT *bc){
  return 0;
}

int bacut_compute_lb(TSP_DATA *d, BACUT *bc){
  int cont = 1, ret;

  while (cont){
    ret = bacut_solve_lp(d, bc);

    if ((ret != OPTIMAL) ||
	(((ret == OPTIMAL) && (!bacut_constraint_generation(d, bc))))) /* Tester ici s'il y a une contrainte a generer */
      cont = 0;

    if (bc->compscount != NULL){free(bc->compscount); bc->compscount = NULL;}
    if (bc->comps != NULL){free(bc->comps); bc->comps = NULL;}
    bc->ncomp = 0; 
  }

  return ret;
}

void bacut_run_rec(TSP_DATA *d, BACUT *bc){
  int i, v;
  double delta, best_delta;
  
  bc->nbnodes++;

  printf("%d/%f\n", bc->nbnodes, (float)bc->best_length); fflush(stdout);

  /* fixer les variables */
  bacut_set_bounds(d, bc);

  if (bacut_compute_lb(d, bc) == OPTIMAL){ /* Le PL a une solution */
    /* Il ne s'agit pas d'une solution entiere  */
   
    /* On cherche une variable sur laquelle brancher */
    v          = -1;
    best_delta = 1.0;
    for (i = 0; (i < d->nbn*(d->nbn-1)/2) /* && (v == -1) */; i++)
      if (((float)bc->x[i] > EPSILON) && ((float)bc->x[i] < 1.0)){
	delta = bc->x[i] - 0.5;
	if ((float)delta < 0.0) delta = -delta; 
	
	if ((float)delta < (float)best_delta){
	  v          = i;
	  best_delta = delta;
	}
      }
    
    /* Inutile dans le branch-and-cut */
    if (v == -1){
      for (i = 0; (i < d->nbn*(d->nbn-1)/2) && (v == -1); i++)
	if (((float)bc->bd[i] == -1))
	  v = i;
    }

    if (v != -1){
      if ((float)bc->x[v] > 0.5){
	bc->bd[v] = 1;
      
	bacut_run_rec(d, bc);
	
	bc->bd[v] = 0;
	
	bacut_run_rec(d, bc);
      }
      else{
	bc->bd[v] = 0;
      
	bacut_run_rec(d, bc);
	
	bc->bd[v] = 1;
	
	bacut_run_rec(d, bc);
      }
      
      bc->bd[v] = -1;
    }
  }
}

void bacut_run(TSP_DATA *d, BACUT *bc){
  int i;
  FILE *f;

  bacut_build_lp(d, bc);
  
  for (i = 0; i < d->nbn*(d->nbn-1)/2; i++)
    bc->bd[i] = -1;
  bc->best_length = MAXDOUBLE;
  bc->nbnodes     = 0;

  bacut_run_rec(d, bc);

  delete_lp(bc->lp);
  
  f = fopen("out.txt", "w");
  fprintf(f, "Meilleure longueur : %7.1f en %d noeuds\n", bc->best_length, bc->nbnodes);
  fclose(f);
}

